st_buffer calculates buffers - polygons encompassing all points up to the specified distance from a given geometry.
We can for example calculate 250km buffers around Avila and Valencia.
To be precise about our distances, we will use EPSG codes
# Transform your data to a suitable projected CRS that uses meters as units# For example, you can use a UTM zone that covers your area of interest# Here, I'm assuming your data is around Avila, Spain, so I'm using UTM zone 30N (EPSG:2062)# See https://epsg.io/ and look for Spain# You should choose an appropriate UTM zone for your data. Note that otherwise, this would be in degreesavila_ctr_sp <-st_transform(avila_ctr, crs =2062)valencia_ctr_sp <-st_transform(valencia_ctr, crs =2062)#Creating the buffer at 250km. Note that unit is in meters: st_crs(avila_ctr_sp) - LENGTHUNIT under Datumavila_250 =st_buffer(avila_ctr_sp, dist =250000)valencia_250 =st_buffer(valencia_ctr_sp, dist =250000)#Going back to the original crs in degreesavila_250 <-st_transform(avila_250, crs =st_crs(esp2))valencia_250 <-st_transform(valencia_250, crs =st_crs(esp2))
Buffers
And this is how we plot the shape file for the buffers that we have just created:
ggplot()+geom_sf(data=esp2)+geom_sf(data=avila_ctr, color ="red", size=3)+geom_sf(data=valencia_ctr, color ="red", size=3)+geom_sf(data=avila_250, color ="red", fill=NA)+geom_sf(data=valencia_250, color ="red", fill=NA)+#geom_sf(data=l, color = "red")+theme_void() +coord_sf(xlim =c(min_lon_x - error, max_lon_x + error), ylim =c(min_lat_y - error, max_lat_y + error))
Buffers
Another type of operation is buffer
st_buffer calculates buffers - polygons encompassing all points up to the specified distance from a given geometry.
We can for example calculate 250km buffers around Avila and Valencia:
Geometry Difference from Pairs
There are four geometry-generating functions that operate on pairs of input geometries:
st_intersection
st_difference
st_sym_difference
st_union
Geometry Difference from Pairs
Geometry Difference from Pairs
Geometry Difference from Pairs
Geometry Difference from Pairs
Geometry Difference from Pairs
Geometry Difference from Pairs: union
We might be interested in the total area that is within 250km from both Anvila and Valencia
We can do this in the following way:
union_area<-st_union(avila_250, valencia_250)ggplot()+geom_sf(data=esp2)+geom_sf(data=avila_ctr, color ="red", size=3)+geom_sf(data=valencia_ctr, color ="red", size=3)+geom_sf(data=union_area, color ="red", fill=NA)+theme_void() +coord_sf(xlim =c(min_lon_x - error, max_lon_x + error), ylim =c(min_lat_y - error, max_lat_y + error))
Geometry Difference from Pairs: union
We might be interested in the total area that is within 250km from both Anvila and Valencia
We can do this in the following way:
Geometry Difference from Pairs: union
Remember that st_union can also be applied to an individual layer,
This returns a dissolved union of all geometries:
Vector layer aggregation
Sometimes, we may not want to dissolve all features into a single geometry
Sometimes, we may want to dissolve by an attribute
Vector layer aggregation
Let us assume that we did not have the polygons for provinces in Spain.
For example, we would be interested in calculating density of railways within provinces
Show the code
#Step1: Reading the polylinesrails <-st_read(dsn="./data/europe-railways-shape/railways.shp", quiet =TRUE)#Step2: Calculating the length of rails in mlengths <-st_length(rails)#Step3: Uniting the linesmerged_line <-st_union(rails)#Step4: Summing up the lengths of the original linestotal_length <-sum(lengths)#Step5: Creating a new sf object with the merged line and total lengthresult_sf <-st_sf(geometry = merged_line, total_length = total_length)ggplot()+geom_sf(data=esp2)+geom_sf(data=result_sf, color ="red")+theme_void() +coord_sf(xlim =c(min_lon_x - error, max_lon_x + error), ylim =c(min_lat_y - error, max_lat_y + error))+ggtitle("Railways in Spain")
Vector layer aggregation
Here is how we calculate the railway density by polygon
#Step1: Performing the intersection between rails and polygonsintersection <-st_intersection(rails, esp2)#Step2: Calculating the length of the intersecting polylinesintersection$length <-st_length(intersection)#Step3: Saving only relevant varsintersection2<-subset(intersection, select =c("NAME_2", "length"))#Step4: Dropping the geometryintersection2<-st_drop_geometry(intersection2)#Step5: Dividing the length by 1000 to obtain kmintersection2$length<-as.numeric(intersection2$length)/1000#Step6: Performing the left join (left = polygons; right = lines)esp2<-left_join(esp2b, intersection2, by =c("NAME_2"="NAME_2"))#Mappingggplot(esp2, aes(fill=length)) +geom_sf() +geom_sf_text(aes(label = NAME_2), colour ="white", size=2)+theme_void() +coord_sf(xlim =c(min_lon_x - error, max_lon_x + error), ylim =c(min_lat_y - error, max_lat_y + error))+scale_fill_viridis_c(option ="viridis")
Vector layer aggregation
Conclusion
We have covered a veriety of sf functions:
perform geometric calculations: areas, changing units
